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Abstract 

Static and dynamic structure factors and various transport coefii- 
cients are computed for a Lennard- Jones model of a binary fluid (A,B) 
with a symmetrical miscibility gap, varying both temperature and rela- 
tive concentration of the mixture. The model is first equilibrated by a 
semi-grandcanonical Monte Carlo method, choosing the temperature and 
chemical potential difference A/i between the two species as the given in- 
dependent variables. Varying for A/i = the temperature and particle 
number A'^ over a wide range, the location of the coexistence curve in the 
thermodynamic limit is estimated. Well-equilibrated configurations from 
these Monte Carlo runs are used as initial states for microcanonical Molec- 
ular Dynamics runs, in order to study the microscopic structure and the 
behavior of transport coefficients as well as dynamic correlation functions 
along the coexistence curve. Dynamic structure factors Sai3{q,t) (and the 
corresponding static functions Sai3{q)) are recorded {a, 13, £ A,B), q being 
the wavenumber and t the time, as well as the mean square displacements 
of the particles (to obtain the self-diffusion constants Da, Db) and trans- 
port coefficients describing collective transport, such as the interdiffusion 
constant and the shear viscosity. The minority species is found to diffuse 
a bit faster than the majority species. Despite the presence of strong 
concentration fluctuations in the system the Stokes-Einstein relation is a 
reasonable approximation. 
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1 Introduction 



While the interplay between static structure and dynamical properties of simple 
fluids is a problem intensively studied both by analytic theory "T','3','3','5','F and 
simulations (5.,.6,, 7, 8..9.if 0). and also extensions to rather complex fluids such as 
polymers are already given [9111)11 1112113114] , somewhat less attention has been 
paid to simple liquid mixtures: while static properties of liquid-liquid miscibility 
have been analyzed in some detail both by analytic theory (see e.g. ^^1) and 
simulations (see e.g. ^MI^M!SMMB^^!^MMM^W}!^M^W!I^2!^ 
|33" 34 '35''36'), dynamic properties of fluid mixtures have been investigated by 
computer simulation only occasionally (e.g. 5 , 35„37..38.^39..40..41.,42j ). In most 
of these studies, the precise phase behavior of the used model is not known, and 
also very small system sizes were used, and hence the proper interpretation of 
the results is difficult. 

The present work aims to contribute towards filling this gap by studying 
both phase behavior and static structure in conjunction with dynamic corre- 
lations and transport phenomena. We address the simplest possible case, a 
symmetric binary Lennard- Jones mixture with a miscibility gap. Extensions of 
our approach to asymmetric Lennard- Jones mixtures and models suitable for 
a realistic description of binary mixtures of liquid metals (such as TO!) will be 
given in future work. We note that a specific asymmetric binary Lennard-Jones 
mixture with 80% A and 20% B particles is a very popular model for the study 
of the dynamics of glass-forming fluids |43II44[I15] . However, the static phase 
diagram of that model system is not yet known. 

The investigation of a symmetric mixture is of particular interest, because 
one can clearly distinguish in this case long-ranged correlations in the vicinity 
of the critical point from the different local ordering of A- and B- particles that 
drives the unmixing transition well below the critical temperature T^: At the 
critical point a mixture of 50% A- and 50% B-particles is formed, i.e. xb = 
1 ~ XA ^ NbUNa + Nb) = 0.5 {Na and Nb being the number of A and B 
particles in the simulation box, respectively), and due to the symmetry of the 
model static A-A and B-B correlations are identical in the latter case. Well 
below Tc at coexistence one of the species, say B, is the minority species and a 
B particle prefers a larger number of neighboring B particles than in the one- 
phase region at the same concentration of B particles. We shall discuss below 
how this interplay of long-ranged correlations in the vicinity of Tc and local 
ordering far below along the coexistence states affects the dynamics of the 
mixture. 

For a symmetric mixture the phase diagram and other static properties 
are most easily computed with Monte Carlo (MC) methods using the semi- 
grandcanonical ensemble rather than Molecular Dynamics (MD) methods. This 
type of MC method has been used with great success to establish binary alloy 
phase diagrams in the solid phase, both for lattice j46ll47ll^l49| and off-lattice 
(50, models. It was also used to study lattice models of polymers [16 25 2 7j^ . 
and, in conjunction with finite size scaling methods, to investigate the crit- 
ical behavior of symmetric binary Lennard-Jones fluids '32'. It is also suit- 
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able for asymmetric models of fluid mixtures [171135) . In this work we combine 
semi-grandcanonical MC simulations with microcanonical MD runs (i.e. at con- 
stant energy E and constant volume V). The former simulations yield well- 
equilibrated configurations at concentrations xb = 1 — a^A and temperatures 
T along the two-phase coexistence curve of the system and thus, these config- 
urations can be used as initial states for the MD to study the structure and 
dynamics at coexistence and also in the vicinity of coexistence by changing the 
temperature or the composition. 

The rest of the Paper is organized as follows. In the next section we shall 
briefly comment on the model and the used simulation techniques. A precise 
definition of the quantities that are computed is given in Sec. 3. Then we present 
the results for the structure (Sec. 4) and dynamics (Sec. 5) of our model that 
have been obtained from MD runs. Finally, Sec. 5 presents our conclusions and 
gives an outlook on further work. 



2 Model and simulation methods 

We consider a binary mixture of point-like particles interacting with Lennard- 
Jones potentials (a,/3 £ A,B), 

(1) 

where {vi} are the positions of the particles, and the Lennard-Jones parameters 
Ea/i are chosen as follows 

CAA = ""BB = CTAB = Cr, ffAA = EBB = £, S — Eab/^- (2) 

We measure all lengths in units of cr = 1, and we also choose Boltzmann's 
constant = 1. The temperature T will be measured in units oi e = 1, and 
(in the MD part) the masses of the two particle species are also chosen the same, 
rUA = "TUB = 1- The potential is truncated and shifted at Vij = 2.5 a. 

The same model has been studied earlier by Wilding [32] (note that in this 
work the potential was also truncated at rij = 2.5 ct but not shifted, and thus 
minor differences in the properties between our model and that of Ref. |32| must 
be expected). In the latter work, also the density of the binary fluid was varied 
over a wide range, < p < 0.7, and it was found for 5 — 0.7, that the A-line 
of critical points for fluid-fluid phase separation ends at a critical end point 
at a density of about pa^ w 0.6 at the coexistence curve of liquid-gas phase 
separation. At sufficiently smaller values of 6 than 0.7 the A line moves towards 
the liquid-gas critical point where it merges into a tricritical point and one may 
expect that this happens at densities below pa^ « 0.6 In the present work, 
we are not interested in the liquid-gas transition of the model at all, and hence 
work at a much higher density, namely p = 1.0. For this density we stay in 
the fluid phase for the whole temperature range of interest, i.e. 1 < T < 1.8. 
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Choosing 6 = 0.5, the critical temperature of unmixing is Tc « 1.638 ± 0.005 at 
this density, see Fig. [21 below. 

The equilibration of the model is done as follows. We start from random 
configurations with an equal number of A and B particles. The particles are 
then displaced by means of a standard Monte Carlo scheme in the NVT ensem- 
ble, i.e. at constant total particle number N, constant volume V, and constant 
temperature T, rejecting or accepting a trial move according to the standard 
Metropolis criterion. Thereby, the length of a trial displacement is a randomly 
chosen number in the interval [— (7/20, (t/20]. After 100000 MC steps the semi- 
grandcanonical MC is switched on: At the end of each displacement step an 
attempted identity switch (B — > A or A — > B) of iV/10 randomly chosen 
particles is made, which is again accepted or rejected according to a standard 
Metropolis criterion, where both the energy change AE and the change of the 
chemical potential ±A/z needs to be taken into account in the Boltzmann fac- 
tor [JliniEIlEniE] ■ Each configuration after this Monte Carlo move (whether 
the change was accepted or not) is counted as a state for the Monte Carlo aver- 
aging in the semi-grandcanonical ensemble, which is defined by the independent 
variables N, V, T, and the chemical potential differences A/i. Note that the 
relative concentration xb = N-b/{Nj^ + Nb) = N-^/N is a fluctuating quantity 
and its average value is an output rather than an input to the simulation. Only 
for the choice A/i = the symmetry of the model between A and B dictates that 
for T > Tc we must have (xb) = 1/2, while for T < Tc we obtain for A/i = 
states at the two-phase coexistence curve. This is illustrated in Fig. where 
the concentration distribution function P{xa) is plotted for three temperatures. 
Note that we have done 5 independent runs, each with a length of 400000 MC 
steps in the semi-grandcanonical ensemble, to calculate the distribution func- 
tions in Fig. ^whereby we started the averaging after 100000 MC steps in each 
run. In principle, for a finite system there is no symmetry breaking possible in 
full thermal equilibrium, and hence one should observe always two peaks, one 
at x™'''^^^^ and the other one at x^''''^^^ = 1 — a;A°''^^^ corresponding to the two 
branches of the two-phase coexistence curve. In practice, these two branches 
are separated at low temperatures by a huge free energy barrier in phase space. 
Therefore, as can be seen in Fig.^ at T = 1.2 and T — 1.4 only the A-rich phase 
is observed, and no transition to the B-rich side of the coexistence curve has 
happened yet at all. Only for T = 1.6, i.e. rather close to (which in Fig. is 
evident from the large width of the peaks) , transitions back and forth between 
the two coexisting phases have occurred (although there were not enough of 
those transitions to make P{xa) strictly symmetric around x^'* = 1/2). 

This "ergodicity breaking" due to finite observation time in the context of 
first-order phase transitions is a well-known phenomenon '51J and does not ham- 
per our analysis at all. From data such as shown in Fig.^we hence obtain, when 
we normalize the part of P(a;A) for xa > 1/2 to J P{xA)dxA — 1, both estimates 
for the location of the coexistence curve and the "concentration susceptibility" 
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^coex(2) ^ J 2,^p(2,^)g^2,^ ^ fcfiTx = iV J x\P{xfC}dxA - [x™'"'''^''^^ 



1/2 



V/2 



(3) 

In order to investigate the dynamics of well-equilibrated states precisely at the 
coexistence curve, we store a number of n = 10 independent configurations at 
each of the 5 independent runs which have a number Nx of A particles such 
that the resulting ratio is as close to possible (note that xa is 

"quantized", because iV^, N are integers, but for large N the discreteness of 
xa, which is l/iV, is smaller than the statistical error with which x™"'^*^^'' can be 
estimated, and therefore this effect does not matter) . With these configurations 
we can realize microcanonical MD runs right at the coexistence curve. Note 
that we have used in this work a standard MD code applying the velocity Verlet 
algorithm with a time step 5t = 0.01 {in units of the time to = {ma^ /ASeY^"^}. 

When we work in the one phase region we just heat up the configurations that 
we have gotten from the semi-grandcanonical Monte-Carlo at coexistence and we 
equilibrate the system by MD at the desired temperature Tf. The temperature 
is then kept constant by coupling the system to a stochastic heat bath, i.e. every 
50 MD steps the velocities are replaced by new ones that are chosen randomly 
from a Boltzmann distribution corresponding to the temperature Tf. 

We have described these simple and almost straightforward considerations [5H 
in such detail here, because the literature is still full of examples where the au- 
thors are not aware of the pitfalls presented to equilibration by the presence 
of slowly relaxing variables. In the present problem, the conservation of the 
concentration (in the ensemble where both TVa and A^b are held fixed) together 
with the presence of a large correlation length of concentration fluctuations 
(which for N ^ oo diverges at the critical point and causes the critical slowing 
down |SS1) may lead to a very slow relaxation. Hydrodynamic slowing down is 
eliminated by the choice of an ensemble where the slow variable is not conserved, 
in our case the semi-grandcanonical ensemble. Critical slowing down cannot be 
avoided (apart from the special case of Ising models on lattices, for which clus- 
ter algorithms have been put forward |55ll5bj that eliminate also critical slowing 
down to a large extent). In the present work, we hence do not attempt to study 
the immediate vicinity of the critical point. 

Fig.|21shows the coexistence curve as estimated from our semi-grandcanonical 
method. One can see clearly that for T < 1.4 finite size effects on the coexis- 
tence curve are quite negligible, if A'' > 400. For T = 1.5, data for A^ = 400 
overestimate the width of the two phase region very slightly, while data for A = 
800, 1600 and 3200 agree within statistical errors. For T = 1.6, however, data 
for both A^ = 400 and A^ = 800 overestimate the width of the coexistence curve 
slightly, while the data for both A^ = 1600 and A^ = 3200 still agree. For this 
reason, we have chosen A^ — 1600 as the standard particle number, for which 
all analysis of static and dynamic correlations have been made (Sees. 3, 4). We 
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have estimated the critical temperature Tc from power law fits of the form 

/(xb) = 1/2 ± = B{1 - T/T,f (4) 

with the critical exponent (3 — 0.33 which corresponds to the universality class 
of the three-dimensional Ising model. B and Tc have been used as adjustable 
parameters. For TV = 400 the fit yields Tc ~ 1.666 ± 0.005, while for all larger 
values of N the result is Tc « 1.638 ± 0.005 (see Fig.E)). We emphasize at this 
point, however, that we deliberately do not address the problem of analyzing the 
critical behavior of transport coefficients in our study: for such a purpose one 
would need to work with system sizes N that are several orders of magnitude 
larger, and extremely accurate data in the temperature range 1.6 < T < 1.7 
would be needed. This interesting problem would require substantially more 
computer resources than were available to us. In addition, we would have to 
locate Tc much more accurately with a finite size scaling analysis [^ 11611^1271 
1^^32115011511 . Even though we do not study critical behavior, our choice of N is 
much larger than most of the choices made in previous work on the dynamics 
of mixtures (e.g. Vogelsang and Hoheisel [35] work with N = 108 and N = 256 
particles, and most data of Asta et al. [3Sj refer to = 500 particles). 



3 Static and dynamic quantities 

We now proceed to define the quantities that are computed in our simulations. 
With respect to static quantities, we have monitored the standard pairwise 
radial distribution functions gap{r) between the diflFerent pairs {a,/3 G (A,B)} 
of particles P[21[ni[l[Sl[ni0[H[ , 



\i=i j=i 



The prime in the second sum means that i — j has to be left out if a = /3. 

In terms of the gaf3{'r), partial structure functions Sapiq) are defined as 
follows [2|: 

OO 

gap{r) Anr'^dr. (6) 



with Xa — Na/N {a e A, B). In order to compute the partial structure factors 
we used the formula 



Safsiq) = J^'^Yl (e^P(*9 ■ ^ki)) , (7) 



fc=i 1=1 

and performed the angular integration numerically to obtain Sa/s^q). 
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It turns out that it is more useful to form the following linear combinations 
of the radial distribution functions ga0 and the partial structure factors Safi 
which describe correlations in number density {gnn and S'nn), concentration 
correlations {gcc and Sec), and cross correlations between number density and 
concentration (gnc and Snc) Pl l57ll58j . The real space functions are given by 

gnn{r) = x%gAA{r) +2xAXBgABiq) + x^gBB{r) , (8) 
gcc{r) = x\xl[gAA{r) + gBB{r) - 2gAB{r)] , (9) 
gncir) = xaxb [XAgAAir) - XBgBBir) + {XB - XA)gAB{r)] , (10) 

and the functions in reciprocal space by 

Snniq) = S'AA(g) + 2S'AB(g) + 5'BB(g) , (11) 
Scciq) = x'^SAAiq) + XASBB{q) -'^XAXBSAsiq) i (12) 
Snciq) = XBSAAiq) - XASBBiq) + {XB - XA)SAB{q) ■ (13) 

We will see below that strong concentration fluctuations are observed if one 
approaches the critical point of unmixing (and this happens not only in the 
immediate vicinity of the critical point). One can try to describe Scciq) at 
small q by the Ornstein-Zernike form, 



5cc(g) = Y^ . (14) 



where the "susceptibility" x is given by Eq. © and ^ has the meaning of a 
static correlation length. 

Turning to dynamic quantities, the most straightforward quantity to consider 
is the mean square displacement of tagged particles 

9A{t) = ([r;,A(0) - nMt)?) , 9B{t) = ([r,.B(0) - rMt)?) , (15) 

where it is understood that the average (• • • ) includes an average over all par- 
ticles of type A or B, respectively. Another quantity which also monitors the 
motion of individual particles are the incoherent intermediate structure func- 
tions F^^\q,t) and F^^\q,t), 

F^^\q,t) = ^(exp{-^g.[rl(0)-r;(^)]}) , (16) 

i GA 

Fl'^Hq,*) - 5](exp{-^g-.[r^,■(0)-r;■(^)]}) , (17) 

where the sum runs over all A particles {Eq. (|16|l } or B particles {Eq. H17(l }. 
respectively. Note that due to the isotropy of the fluid, Eqs. (I16|l . H17|l cannot 
depend on the direction of the scattering vector q. When the displacements of 
the particles are Gaussian distributed, Eqs. (|16|l . (|17|l reduce to 



F^^\q,t) 



exp 



^q^gAit) 
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Fl^'Kq^t) 



exp 



Iq^gBit) 



(18) 
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Noting further the Einstein relations 



(19) 



where Da and Db are self-difFusion constants of the particles, we expect that 
the asymptotic decay of Fi^-* (g, t) and f]^'' {q, t) can be described by a simple 
exponential variation with time, 



FiA)(g,t) cxexp[-</rA((?)] 
and Eqs. (23-1101 then imply 

lim[TA(g)g^]"^ = -Da, 



FiB)('?,i) «exp[-VTB(<z)] 



lim[rB(g)q^ 

9^0 



(20) 



(21) 



In addition, we introduce partial coherent intermediate structure functions 
Safjil^i) defined as [a,/3 G (A,B)] 



(22) 



1=1 j=i 
iea j e IS 



and from these functions we can define number density and concentration cor- 
relation as well as cross correlation functions {5„„(g,i), Scci,q,t) and 5„c(<Z,i)} 
by forming exactly the same type of linear combinations as written in Eqs. 
ifOjl . but with the Sajsiq, t) rather than their static counterpart Sa/siq)- In order 
to compare the time correlation functions Sapiq, t) for different values of q it is 
convenient to consider the normalized functions 



F^q) 



Safj{q,t) 



(23) 



The normalized functions Fnn{q,t), Fcc{q,t), and Fnc{q,t) are defined in a sim- 
ilar way by normalizing S'„„(q,t), Scc{q,t) and Snc{q,t) by the corresponding 
static correlation function. 

Further quantities of interest are the instantaneous values of the components 
of the pressure tensor (x, y, z denote the Cartesian components, Vi is the velocity 
of the i'th particle) 



N 



i=l 



■'^jl-Fyilfi Tj 



(24) 



F being the force acting between particles From the autocorrelation func- 
tion of (Jxyit) one gets the shear viscosity ry of the fluid by a Green-Kubo for- 
mula PE], 



1 

VkeT 



dt {(Txy{0)axy{t)) . 



(25) 
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Finally we mention the Green-Kubo formula for the interdiffusion constant 
(note that we have chosen the masses of all particles equal to each other here) |2] 



oo 



NScciO) 



1 



/ 



(26) 







where the current J' 



(t) is defined as follows: 




(27) 



1=1 1=1 



One can easily calculate I?int for an ideal mixture [such a mixture is formed 
if the particles are labeled but are otherwise identical and thus the entropy of 
mixing is equal to —NksixA^ogxx + 2;BlogXB)]. In this case cross correlations 
in the autocorrelation function for J^"' in Eq. (|26|1 vanish and the interdiffusion 
constant can be expressed by the self-diffusion constants 2 , 



It has been shown by MD simulations that Eq. H28|) is a good approximation for 
mixtures of Lennard- Jones fluids in the one- phase region ,38n59i| . However, for 
ideal mixtures on a rigid lattice (where A, B-atoms hop with jump rates cta, cb 
to vacant sites, assuming that a small number of vacancies is present) Eq. I|28|l 
does not hold |60j . 

4 Static properties of the symmetrical binary 
fluid mixture 

In this section we analyze the structure of the symmetric Lennard- Jones mix- 
ture. Results are shown along the coexistence line and along the different paths 
in the one-phase region that are indicated in Fig. |2 

In Fig. 131 we present the pair correlation functions gapif) for the three states 
at temperatures T = 1.2, 1.4, and 1.6 at the coexistence curve (cf. Fig.lJJ. These 
data show a typical normal fluid behavior in all cases, as expected. Although it 
would be hard to recognize from these partial radial distribution functions that 
one approaches the critical point of unmixing in a binary fluid, the comparison 
of the different correlations shows nontrivial features. At each temperature, 
the first peak in gAB{f) is at a slightly smaller distance than in gAAi''') and 
gBBif)- Moreover, the second peak in 5bb(^) is shifted towards smaller distances 
compared to the other two functions and this effect is more pronounced the lower 
the temperature is, i.e. the farther one is away from the critical point. The 
explanation of this feature is simple: The B particles are the minority species 
at the considered state points and thus, at small concentrations xb it is very 
likely that one finds an A particle between two next-nearest B particles forming 



Ant = xbDa + xaDb ■ 



(28) 
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a B-A-B sequence. Since the distance between nearest A-B neighbors is smaller 
than that between nearest A-A and B-B neighbors one may expect that a large 
amount of such B-A-B sequences, where the distance between next-nearest B 
particles is smaller than in a B-B-B sequence, leads to the shift in the second 
peak of gBB{r)- Certainly, the difference between gAA and 5bb becomes less 
pronounced the closer one approaches Xb = 0.5 (and thus the critical point) 
where, due to the symmetry of our model, gAA and ijbb are identical. 

The temperature dependence of gAAif) and gBBif) in the one-phase region is 
shown in Fig.^jfor the constant concentration xb = 0.10375 in the temperature 
range 1.4 < T < 1.8 (T = 1.4 corresponds to a state on the coexistence curve). 
We see that the temperature dependence is very weak for gAAir) whereas we 
observe a significant increase of the amplitude of the first peak in gBsir) de- 
creasing the temperature towards the coexistence line. The latter effect is even 
more pronounced in the coordination number function Zaair) that is shown in 
the insets of Fig. 0] This function is defined by 



which gives the number of particles of type a surrounding a particle of type a 
within a distance r' < r. In the considered temperature range 2:bb('') changes 
at r = 1.4 (i.e. around the location of the first minimum) from zbb = 1-37 for 
T = 1.8 to zbb = 1.55 for T — 1.4. By a closer inspection of ZAA{r) one finds 
that the changes in this quantity are of the same order. Thus there are only 
minor changes in the local order if one approaches the coexistence line from 
above. 

We now proceed to discuss the behavior of the structure factors 5„„(g), 
Scc{q), and Snciq) {Eqs. C^-CSl} along the coexistence line and for the latter 
two quantities also the corresponding functions in real space {Eqs. l(^- H10|l l. see 
Fig. 13 The number density structure factor <S'n„(g) behaves exactly as expected 
for any dense simple liquid, and in particular Snnil 0) is very small, as 
expected for liquids that are almost incompressible. The height of the first peak 
is around 3, as is typical for fluids at temperatures somewhat higher than the 
melting temperature, and there is a pronounced second peak, due to the rather 
regular close packing of atoms in a dense fluid. 

In contrast, Scciq) at T = 1.2 is very small and structureless throughout: in 
symmetrical binary mixtures, there is very little coupling between density and 
concentration fluctuations. For T = 1.4, we see already some enhancement of 
Sccil) at small q, reflecting the growing concentration fluctuations. For T = 1.6, 
at small q a rather dramatic growth of Scc{q — 0) signals the proximity of the 
critical point. For small q we have fltted the data by the Ornstein-Zcrnikc form 
Eq. H14|l (bold lines in Fig.|SlD), where we have used the correlation length ^ as 
an adjustable parameter and we have estimated the susceptibility x by means of 
Eq. Q from the distribution function P{xa) (the values of S'cc(O) = ksTx are 
shown as crosses in Fig. (SJ)) . From the fits we recognize that the data for small 
q can be described by Eq. H14I) within the statistical errors (note that there is 
only one adjustable parameter, the correlation length ^). From the fits we find 




a E [A, B] 



(29) 
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a growing correlation length ^ with increasing temperature towards T^. This 
behavior can be also inferred from gcc{r) (see inset of Fig.|S|D) at T = 1.6 which 
has a rather pronounced tail for r > 2. Also remarkable in gcc{f) is the negative 
peak around r k, 0.9 which stems from the fact that these distances are avoided 
by nearest A-A and B-B neighbors but are typical for nearest A-B neighbors. 

Particularly interesting is the structure factor Sndl)- In a mixture with 
A/Lt = at T > Tc, where then xa = a^B = 1/2, we expect that Snc{<l) = due 
to the perfect symmetry between A and B. Along the coexistence curve, the 
spontaneous symmetry breaking between A and B invalidates this argument, 
and hence a nonzero Snc{q) is possible. As we see in Fig. this is obviously 
the case. The origin of the oscillations in Sndq) can be most easily understood 
by means of its counterpart in real space, gnc{r) (see inset of Fig. |SJ;). The 
most pronounced feature in the latter quantity is a peak with negative intensity 
around ri w 0.9. This feature is again due to the fact that the first peak in 
AB (r) is slightly shifted to smaller distances compared to the first peak in ^aa 
and in ^bb such that the intensity of the latter two functions in the interval 
0.85 < r < 0.9 is essentially zero whereas gAB is rapidly increasing for the 
latter distances. This leads to the occurrence of negative correlations in gnc 
and thus to the peak around ri. We emphasize that the scale for 5'„c('z) and 
g„c(r) is clearly very much smaller than for the other structure factors and pair 
correlation functions, respectively. 

Similar data have been taken for two other paths in the phase diagram of 
Fig. 121 namely one at constant temperature, T = 1.5, varying A/i (and hence 
Xb), and one at constant xb = 0.10375. It turns out that the radial distribution 
functions gafsir) always look similar to those that are shown in Fig.|21 and hence 
are not shown here. The same statement applies to <5'„„(q), which exhibits only 
a weak dependence on either temperature or composition. Here we hence show 
only those quantities which have a more pronounced dependence on our control 
parameters, namely Scc{q) and Sndq), Figs. 0. From Fig.|BJi we see that 
the enhancement of Scdl — * 0) persists to rather high temperatures, and also 
Sncil) is only weakly temperature-dependent in the temperature region shown. 
In contrast, when we change the composition at T = 1.5, we find that Sccil) 
quickly loses all structure when xb ^ (Fig. [7^), and a related trend is observed 
in Snc{q)- 

5 Dynamic properties of the symmetrical binary 
fluid mixture 

The detailed description of the structure of our Lennard-Jones model that we 
have done so far is necessary to understand its dynamic properties. The next 
two sections are devoted to the results for time-dependent correlation functions 
and transport coefficients. In Sec. 5.1 we analyze the dynamic properties of 
our system along the fluid-fluid coexistence line whereas in Sec. 5.2 dynamic 
quantities are discussed that have been obtained along the path at xb = 0.10375 
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which starts in the one-phase region and ends at the coexistence Une at T = 1.4 
(see Fig.El). 

5.1 The dynamics at phase coexistence 

We start by showing the incoherent intermediate structure function for the A 
particles, Fs^\q,t) {Eq. ((TCll T. see Fig. |S1 (data for B-particles look very similar 
and therefore are not shown). Note that for q « 7, where the peaks of the 
static structure factors Sapiq) and Snn{q) occur, F"{q,t), a £ (A,B), decays 
on the time scale of about t = 20 and this decay is certainly even faster for 
q = 10. Moreover, the temperature dependence of the F"{q,t) is weak in the 
temperature range under consideration. Such a short structural relaxation time 
as well as a weak temperature dependence indeed are expected for any ordinary 
fluid, which has no tendency to glass formation |4HII44lEKj . 

On the other hand, for small q the decay of Fs'^\q,t), Fs^\q,t) is much 
slower, and this slowing down simply reflects the diffusive behavior predicted in 
Eqs. H18 |l - H21(l . Fig. El shows that this interpretation in fact is nicely consistent 
with the data in that, in agreement with Eq. (|3U|I . the product r^g^ approaches 
the inverse diffusion constant Da for q (the values for D^^ are indicated 
as arrows in Fig. |5Jl. One can further see that the B particles diffuse a bit 
faster than the A particles. Of course, we expect that this dynamic asymmetry 
must vanish when we are at the critical point, and indeed, as one can infer from 
Fig. El the ratio Db/Da decreases to one if the critical point is approached from 
below along the coexistence line. The difference in the diffusion constants has 
a simple reason: In the A rich phase the B particles find mainly A particles as 
nearest neighbors, and since the interaction between particles of different species 
is weaker than for equal particles it is easier for the B particles to escape from 
an A rich neighborhood than from a B rich one. Thus, the effective activation 
energy for the diffusion of B particles is smaller in an A rich than in an B 
rich environment. Note that the behavior of r^g^ at finite q is very different 
from liquids near a glass transition: In such liquids one would see a peak in 
TQ,q^ around a q value that corresponds to the location of the structure factor 
maximum (see Ref. |(jl) and references therein). 

Particularly interesting is the behavior of the collective structure functions 
fFigs. [TUl lll|) . In FAA{q,t) (Fis. I10() we observe again a very fast decay for q 
values around the structure factor maximum q = 7 and at q = 0.992, i.e. at small 
q, the decay occurs on a two orders of magnitude larger time scale. For the latter 
value of q one can clearly recognize that several relaxation processes contribute, 
since the decay emerges in two steps. Note that Fab (<?, t) and Fbb {<!, t) exhibit a 
similar behavior and so are not shown here. In contrast to the partial structure 
functions, in the concentration-concentration correlation function Fcc{q,t) only 
the slow decay is found at g = 0.992 fFig. lll(l . Surprisingly, this is also the case 
at r = 1.2 although, at this temperature, there is almost no structure in the 
corresponding static function Sccil)- As can be also inferred from Fig. 1111 the 
function Fn„(q, t) shows a very fast decay to zero at g = 0.992 (accompanied by 
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oscillations due to acoustic modes) , whereas it has a slightly slower relaxation to 
zero than Fcc{q, t) at q = 7, i.e. the location of the structure factor maximum. 

The relaxation times for the collective correlation functions are of the same 
order as those detected from the self-correlation functions: disappointingly, 
there is no clear indication of critical slowing down when we approach Tc from 
below, since in the temperature region studied the overall Arrhenius-like increase 
of the relaxation time as T is lowered dominates (the Arrhenius-like behavior is 
demonstrated below in the case of the diffusion constants). 



5.2 The dynamics from the one phase region to coexis- 
tence 

The qualitative character of the relaxation functions Fs''^\q,t), Fs^\q,t) and 
FAA{q, t), Fcc{q, t) in the one phase region is rather similar to the behavior found 
along the coexistence curve. Therefore, we do not show these functions in any 
detail here, but proceed immediately to a counterpart of Fig.lHlin Fig.^J where 
now the relaxation times Tx(ci)q^ and 'TB(q)q^ extracted from the incoherent 
intermediate scattering functions at constant xb = 0.10375 are shown as a 
function of wavcnumber q. Again one finds that the minority species diffuses 
a bit faster than the majority. The reason for the larger discrepancies of 
and Ta{q)q^ for g ^ in the case of the B particles is of course due to the fact 
that the statistics is worse for the B particles which is the minority species. 

This difference in the diffusion constants is emphasized in the Arrhenius plot, 
Fig. 1131 where the logarithm of Da is plotted vs. inverse temperature. It is seen 
that a simple thermally activated behavior (which would show up as straight 
lines in this plot) holds only over very restricted temperature regimes which 
might be related to the fact that there are small but significant changes in the 
local structure when approaching the coexistence line (see especially Fig. EJ. 
Also plotted in Fig. [H are the inverse shear viscosity rj ^ and the interdiffu- 
sion constant i'int. For an accurate computation of these quantities we have 
undertaken a much larger effort than for the single particle quantitities: At 
each temperature we averaged over 50 independent runs of 400000 MD steps to 
calculate rj and i^int by means of Eq. H25|l and Eq. (|26|l , respectively. We show 
also in Fig. E| the interdiffusion constant according to Eq. (|28|l that is valid 
for an ideal mixture. Already for T < 2.5 I?int decouples from the ideal mix- 
ture approximation and thus from the self-diffusion constants, and Z?int becomes 
significantly smaller than the Da if one approaches the coexistence curve. 

It is now interesting to check for the Stokes-Einstein relation, which relates 
to the self-diffusion constant of a diffusing spherical particle of diameter ci in a 
fluid of viscosity 77 as 

D = M_. (30) 

The factor 2 in the denominator of Eq. H30|) corresponds to the assumption 
of slip boundary conditions on the surface of the diffusing particle (for stick 
boundary conditions a factor 3 instead of 2 appears) [^21 • 
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Using the data for 77, Dx, Db and invoking Eq. (|30|) one obtains the corre- 
sponding effective Stokes- Einstein diameters c^a , shown in Fig. The data 
for da show that despite the strong concentration fluctuations in the system the 
Stokes-Einstein relation is a good approximation. The finding that the values 
for dA and are different is reasonable since the distance for nearest A-B neigh- 
bors is slightly smaller than that for A-A or B-B neighbors and the B particles 
as the minority species are mostly surrounded by A particles (see Fig.^J; thus, 
the effective hydrodynamic diameter for the B particles is slightly smaller than 
that for the A particles. Note that with the assumption of slip boundary con- 
ditions in Eq. Ij^Ufl we obtain reasonable values around 1 for dA and de whereas 
for stick boundary conditions these values are rather small (around 0.6). 

6 Concluding remarks 

The present paper considers as a simple model system for a binary fluid mixture 
with a miscibility gap a dense Lennard-Jones mixture with fully symmetrical 
interactions, ctaa = cbb = o^ab = o", eaa = £bb = s = 1, S = eab—s = 0.5, and 
also the masses of the two types of particles are chosen the same, = mB = m. 
The aim of the present work was to carry out a feasibility study, where by 
combination of the semi-grandcanonical Monte Carlo technique with Molecular 
Dynamics simulations both static properties, including the phase diagram, and 
dynamic correlation functions and associated transport coefficients are obtained 
simultaneously. It is shown that outside of the critical region of the mixture 
rather small system sizes (such as a total number of iV = 1600 particles) already 
suffice to obtain quantitatively reliable results, while in the critical region (which 
could not yet be studied here) much larger sizes clearly are mandatory. Such 
large sizes would also be needed for the very interesting problem of spinodal 
decomposition of the binary fluid mixture after a quench from the one-phase 
region into the unstable region of the phase diagram (Fig.|2Jl. Since our study is 
the first study of a fluid mixture where both the phase diagram and the relevant 
transport coefficients have been determined simultaneously, such an extension 
to far from equilibrium dynamics clearly would be very interesting. 

Our studies of transport coefficients are distinct from earlier work on similar 
models by the fact the we know precisely where in the phase diagram the con- 
sidered state points are for which transport coefficients have been determined. 
This fact clearly facilitates the proper interpretation of the results. A particular 
interesting finding concerns the asymmetric composition in our otherwise fully 
symmetric model. 

In further work we also plan to study systematically Onsager coefficients 
to test the concepts of phenomenological irreversible thermodynamics for the 
present model and to check "mixing rules" for the interdiffusion coefficient. 
Moreover, we will present extensions to more realistic models of real materials 
which always exhibit asymmetry in the interactions, to allow also a comparison 
to suitable experiments. 
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Figure 1: Distribution function P{xa) for the relative concentration a;A = 
Na/N of A particles for A'' = 1600 particles at density p = 1 and the three 
indicated temperatures. 
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Figure 2: Phase diagram of the symmetrical binary Lennard- Jones mixture in 
the plane of variables temperature T and concentration xb = N-b/N, for density 
p = 1.0 and four choices of N, as indicated. Crosses and plus symbols indicate 
paths where the structure and dynamics is studied at constant concentration 
xb — 0.10375 as a function of temperature and at constant temperature T — 1.5 
as a function of concentration xb, respectively (see Sees. 4 and 5.2). Squares 
with arrows show states at the coexistence curve, which will be analyzed in 
detail in Sees. 4 and 5.1. The full curves show power law fits with Eq. 10}. For 
N — 400 (thin curve) one obtains T^. « 1.666 ± 0.005, while for all larger values 
of N the result is Tc « 1.638 ± 0.005. The broken curves are guides to the eye 
only. The crosses at xs — 1/2 mark the estimates for Tc that result for the 
various choices of N. 
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Figure 3: Radial distribution functions gAAi^) (solid lines), 5ab('') (dashed 
lines), and gBBir) (dashed-dotted lines) along the A-rich part of the coexistence 
curve (Fig.EJ plotted vs. r at the temperatures a) T = 1.2, h) T = 1.4, and c) 
T= 1.6. 
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Figure 4: Radial distribution functions gAAir), part a), and gBsir), part b), 
for the temperatures T = 1.4, 1.5, 1.6, 1.7, and 1.8 at the concentration xb = 
0.10375. The insets show the corresponding coordination numbers Zaa ol e 
[A, B] as a function of r. 
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Figure 5: Different structure factors plotted vs. q for the three indicated tem- 
peratures along the coexistence curve, a) Snn{q), b) Scc{q) (the bold solid hues 
are fits with Eq. El the crosses are estimated of Scc{q — 0) as obtained from 
Eqs. jnj and (|14(l ^. and c) Snc{q)- Note that a logarithmic scale is used on the 
ordinate of part b). The insets in b) and c) show the corresponding functions 
in real space. 
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Figure 6: Structure factors Scciq), part a), and Snc{q), part b), plotted vs. q, 
for the temperatures T = 1.4, 1.5, 1.6, 1.7, and 1.8 (from top to bottom), at 
the concentration xb = 0.10375. For Sec again a logarithmic scale is used. The 
curves are overshifted relative to each other by an amount A = 0.005 (T = 1.8 
is the unshifted curve). The crosses in part a) are estimates of Scc{q = 0) 
as in Fig. \^ and the bold solid lines are fits with Eq. E| The insets show 
the corresponding functions in real space (also overshifted to each other by an 
amount A = 0.001 in part a) and A = 0.005 in part b), again T = 1.8 is the 
unshifted temperature). 
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Figure 7: Same as Fig. but at constant temperature T = 1.5 and three 
concentrations as shown (note that aU the curves are unshifted). 
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Figure 8: Incoherent intermediate scattering function of A particles plotted 
versus time (note the logarithmic scale of time) for the three temperatures 
T = 1.2 (solid lines), T = 1.4 (dashed lines), and T = 1.6 (dashed-dotted lines) 
along the coexistence curve for three different values of q, namely: q = 0.992, 
7, 10 (from right to left). 



27 




Figure 9: Plot of scaled relaxation time T{q)q^ vs. q, for the A-particles (open 
symbols) and the B-particles (filled symbols) for three temperatures along the 
A-rich branch of the coexistence curve. The arrows indicate the value of the 
corresponding inverse self-diffusion constants. 
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Figure 10: Coherent intermediate scattering function FAAiq,t) of A-particles 
plotted vs. time (note logarithmic scale of time) for three temperatures and 
concentrations given from the A-rich part of the coexistence curve. Three tem- 
peratures are shown, T = 1.2, part a), T = 1.4, part b), and T = 1.6, part c). 
Three wavenumbers q are included as indicated. 
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Figure 11: Same as Fig. E|but for Fcc{q,t) (solid lines) and F„„(g, (dashed 
lines) . 
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Figure 12: Scaled relaxation times Tj^{q)q^ for the A-particles (open symbols) 
and T^{q)q^ for the B-particlcs (filled symbols) plotted vs. wavemimbcr q. Three 
temperatures at constant xb = 0.10375 are shown. The arrows indicate the 
values of the corresponding inverse diffusion constants, as extracted from the 
mean square displacements via the Einstein relation. 
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Figure 13: Logarithm of the self-difFusion constants (fiUed symbols) and the 
interdiffusion constant (open diamonds) as well as of the inverse viscosity (open 
triangles) plotted vs. 1/T for a wide range of temperatures for the mixture at 
constant concentration xb — 0.10375. Also plotted is Dint as calculated from 
Eq. H28|l (open circles). The viscosity is multiplied by a factor of 3. The lines 
are guides to the eye. 
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Figure 14: Stokes-Einstein diameters cIa, ds of the diffusing particles plotted 
vs. temperature, for the mixture at constant concentration xb = 0.10375. 
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